﻿using dnAnalytics.LinearAlgebra;
using dnAnalytics.LinearAlgebra.Decomposition;
using dnAnalytics.LinearAlgebra.Solvers.Direct;
using NUnit.Framework;

namespace dnAnalytics.Tests.LinearAlgebra.Solvers.Direct
{
    public abstract class QRSolverTests
    {
        private const string squareMatrixFile = "./TestData/Matrices/random_real_general_array_100.mtx";
        private const string tallMatrixFile = "./TestData/Matrices/random_real_general_array_20_10.mtx";

        protected abstract Matrix GetMatrix(string file);
        protected abstract Matrix GetMatrix(int rows, int columns);
        protected abstract Matrix GetIdentityMatrix(int size);
        protected abstract Vector GetVector(int count);

        [Test]
        public void HouseholderSolveSquareMatrix()
        {
            Matrix squareMatrix = GetMatrix(squareMatrixFile);
            Matrix identity = GetIdentityMatrix(squareMatrix.Rows);
            for (int i = 0; i < identity.Rows; i++)
            {
                identity[i, i] = 1.0;
            }

            QRSolver<Householder> solver = new QRSolver<Householder>();
            Matrix result = solver.Solve(squareMatrix, identity);

            result = squareMatrix*result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.AcceptableError);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        public void HouseholderSolveTallMatrix()
        {
            Matrix tallMatrix = GetMatrix(tallMatrixFile);
            Matrix matrix = GetMatrix(tallMatrix.Rows, 5);
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    matrix[i, j] = (i + 1.0)/(j + 2.0);
                }
            }

            QRSolver<Householder> solver = new QRSolver<Householder>();
            Matrix result = solver.Solve(tallMatrix, matrix);
            Matrix test = (tallMatrix.Transpose()*tallMatrix).Inverse()*tallMatrix.Transpose()*matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(test[i, j], result[i, j], Constants.AcceptableError);
                }
            }
        }


        [Test]
        public void HouseholderSolveSquareVector()
        {
            Matrix squareMatrix = GetMatrix(squareMatrixFile);
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0)/(i + 2.0);
            }

            QRSolver<Householder> solver = new QRSolver<Householder>();
            Vector result = solver.Solve(squareMatrix, vector);
            result = squareMatrix*result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0)/(i + 2.0), result[i], Constants.AcceptableError);
            }
        }

        [Test]
        public void HouseholderSolveTallVector()
        {
            Matrix tallMatrix = GetMatrix(tallMatrixFile);
            Vector vector = GetVector(tallMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0)/(i + 2.0);
            }

            QRSolver<Householder> solver = new QRSolver<Householder>();
            Vector result = solver.Solve(tallMatrix, vector);

            Vector test = (tallMatrix.Transpose()*tallMatrix).Inverse()*tallMatrix.Transpose()*vector;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual(test[i], result[i], Constants.AcceptableError);
            }
        }

        [Test]
        public void GramSchmidtSquareMatrix()
        {
            Matrix squareMatrix = GetMatrix(squareMatrixFile);
            Matrix identity = GetIdentityMatrix(squareMatrix.Rows);
            for (int i = 0; i < identity.Rows; i++)
            {
                identity[i, i] = 1.0;
            }

            QRSolver<GramSchmidt> solver = new QRSolver<GramSchmidt>();
            Matrix result = solver.Solve(squareMatrix, identity);

            result = squareMatrix * result;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    if (i == j)
                    {
                        Assert.AreEqual(1, result[i, j], Constants.AcceptableError);
                    }
                    else
                    {
                        Assert.AreEqual(0, result[i, j], Constants.AcceptableError);
                    }
                }
            }
        }

        [Test]
        public void GramSchmidtTallMatrix()
        {
            Matrix tallMatrix = GetMatrix(tallMatrixFile);
            Matrix matrix = GetMatrix(tallMatrix.Rows, 5);
            for (int i = 0; i < matrix.Rows; i++)
            {
                for (int j = 0; j < matrix.Columns; j++)
                {
                    matrix[i, j] = (i + 1.0) / (j + 2.0);
                }
            }

            QRSolver<GramSchmidt> solver = new QRSolver<GramSchmidt>();
            Matrix result = solver.Solve(tallMatrix, matrix);
            Matrix test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * matrix;

            for (int i = 0; i < result.Rows; i++)
            {
                for (int j = 0; j < result.Columns; j++)
                {
                    Assert.AreEqual(test[i, j], result[i, j], Constants.AcceptableError);
                }
            }
        }


        [Test]
        public void GramSchmidtSquareVector()
        {
            Matrix squareMatrix = GetMatrix(squareMatrixFile);
            Vector vector = GetVector(squareMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }

            QRSolver<GramSchmidt> solver = new QRSolver<GramSchmidt>();
            Vector result = solver.Solve(squareMatrix, vector);
            result = squareMatrix * result;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual((i + 1.0) / (i + 2.0), result[i], Constants.AcceptableError);
            }
        }

        [Test]
        public void GramSchmidtTallVector()
        {
            Matrix tallMatrix = GetMatrix(tallMatrixFile);
            Vector vector = GetVector(tallMatrix.Rows);
            for (int i = 0; i < vector.Count; i++)
            {
                vector[i] = (i + 1.0) / (i + 2.0);
            }

            QRSolver<GramSchmidt> solver = new QRSolver<GramSchmidt>();
            Vector result = solver.Solve(tallMatrix, vector);

            Vector test = (tallMatrix.Transpose() * tallMatrix).Inverse() * tallMatrix.Transpose() * vector;

            for (int i = 0; i < result.Count; i++)
            {
                Assert.AreEqual(test[i], result[i], Constants.AcceptableError);
            }
        }
    }
}